% 创建两个关于时间t变化的信号
syms t S1 S2;
%% 有用信号
S1t = exp(-abs(t));  % 双边指数函数
F1 = fourier(S1t);   % 对信号S1t进行傅里叶变换到频域(w)上去

% Sa1的时域波形和频谱图
figure(1);
fplot(S1t);   % S1t就是fplot要画的函数
axis([-5,5,-0.5,1.2]);
title('S_a_1(t)的时域波形');
xlabel('t');
grid on;     % 显示网格线
figure(2);
fplot(F1);
title('S_a_1(t)频谱波形');
axis([-5,5,0,2.5]);
xlabel('w');
grid on;

%% 干扰信号
S2t = 10*exp(-abs(t)).*cos(10*t);    % 20dB就是10倍，并将频谱进行搬移
F2 = fourier(S2t);                   % 使有用信号和干扰信号的频谱不重叠
% Sa2t的时域波形和频谱图
figure(3)
fplot(S2t);
axis([-3,3,-11,11]);
title('S_a_2(t)时域波形');
xlabel('t');
grid on;
figure(4);
fplot(F2);
title('S_a_2(t)频域波形');
axis([-20,0,0,15]);
xlabel('w');
grid on;

%% 合成信号，将两个信号叠加
Xat = S1t + S2t;
F3 = fourier(Xat);
% Xat的时域波形和频谱图
figure(5);
fplot(Xat);
axis([-2.5,2.5,-11,11]);
title('X_a(t)的时域波形');
xlabel('t');
grid on;
figure(6);
fplot(F3);
title('X_a(t)频域波形');
axis([-30,0,0,12]);
xlabel('w');
grid on;

%% 对上面三个时域信号采样
N = 100;       % 100个采样点
n = 0:1:N-1;
fs = 30;   % 采样频率为30Hz
Ts = 1/fs;
S1n = exp(-abs(n*Ts));    % t = n*Ts 第n+1个采样点，采样点之间间隔为Ts
S2n = 10*exp(-abs(n*Ts)).*cos(10*n*Ts);
Xan = S1n+S2n;
S1k = abs(fft(S1n,N));    % 快速傅里叶变换
S2k = abs(fft(S2n,N));    % 采样后得到的信号转换到频域(k)上
Xak = abs(fft(Xan,N));
% 画上述三个采样后的离散信号的时域和频域波形
% 画出S_1(n)的时域和频域上的图形
figure(7)
stem(n,S1n,'.');
grid on;
title('S_1(n)的时域波形');
xlabel('n');
ylabel('S_1(n)');
figure(8)
stem(n,S1k,'.');
title('S_1(n)的频谱图');
xlabel('k');
ylabel('S_1(k)');
grid on;
% 画出S_2(n)的时域和频域上的图形
figure(9)
stem(n,S2n,'.');
grid on;
title('S_2(n)的时域波形');
xlabel('n');
ylabel('S_2(n)');
figure(10)
stem(n,S2k,'.');
title('S_2(n)频谱图');
xlabel('k');
ylabel('S_2(k)');
grid on;
% 画出X_a(n)的时域和频域上的图形
figure(11)
stem(n,Xan,'.');
grid on;
title('X_a(n)的时域波形');
xlabel('n');
ylabel('Xa(n)');
figure(12)
stem(n,Xak,'.');
title('X_a(n)的频谱图X_a(k)');   % X_a(n)的频谱图
grid on;

%% 设计巴特沃斯滤波器
wp = 3;
ws = 7;
Wp = (2/Ts)*tan(wp/2);   % 设置归一化通带和阻带截止频率
Ws = (2/Ts)*tan(ws/2);
Rp = 1;
Rs = 40;  % 阻带最小衰减为40dB
% 双线性变换法
[N,Wc] = buttord(Wp,Ws,Rp,Rs,'s');   % 计算模拟滤波器的阶和3dB截止频率
[B,A] = butter(N,Wc,'low','s');      % 计算模拟滤波器系统函数分子、分母多项式系数
B
A
W = 0:0.001:50;
Hs = freqs(B,A,W); % 计算幅频响应
[Bz,Az] = bilinear(B,A,fs);      % 用双线性变换法法将模拟滤波器转换成数字滤波器
Hz = freqz(Bz,Az,W);     % 返回频率响应
phi1 = angle(Hs);
phi2 = angle(Hz);
Hs = 20*log10(abs(Hs));      % 幅度，单位为dB
% 画出滤波器的幅频响应和相频响应曲线
figure(13);
plot(W,Hs);
title('模拟滤波器的幅频响应');
xlabel('w');
ylabel('|Hw|');
grid on;
figure(14);
plot(W,phi1);
title('模拟滤波器的相频响应');
xlabel('w');
ylabel('angle(Hw)');
grid on;

% 画出数字滤波器的幅频响应和相频响应曲线
figure(15);
plot(W/pi,abs(Hz));
title('数字滤波器的幅频响应');
xlabel('w');
ylabel('|Hw|');
grid on;
figure(16);
plot(W,phi2);
title('数字滤波器的相频响应');
xlabel('w');
ylabel('angle(Hw)');
grid on;



%% 用上面设计好的巴特沃斯滤波器滤波
yn = filter(Bz,Az,Xan);     
n = 0:length(yn)-1;
Yk = abs(fft(yn,length(yn)));
% 画滤波之后的时域和频域波形
figure(17)
plot(n,yn);
title('滤波之后的y_n时域波形');
xlabel('n');
ylabel('yn');
grid on;
figure(18)
stem(n,Yk,'.');
title('滤波之后的y_n的频谱');
xlabel('w');
ylabel('Yk');
grid on;

 
